Add standalone reference tests for cellular ODE models - #627
Conversation
|
The Regazzoni reference generator currently uses the authors’ C++ implementation at a pinned commit, so regenerating it requires a local clone of I did consider making it more self-contained, either by optionally fetching the pinned source automatically or by writing an independent Python version similar to the TP06 generator. I also wasn’t sure whether bundling or automatically fetching an external reference implementation is something we want to do in svMultiPhysics, so I didn’t want to make that choice without asking. I’d be happy to change the workflow if you think a more self-contained approach would be preferable. |
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #627 +/- ##
==========================================
+ Coverage 72.86% 73.08% +0.22%
==========================================
Files 258 266 +8
Lines 39500 39980 +480
Branches 6730 6797 +67
==========================================
+ Hits 28780 29221 +441
- Misses 10477 10516 +39
Partials 243 243 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
michelebucelli
left a comment
There was a problem hiding this comment.
Thank you @samibismar! This is a very useful contribution. I have left some comments below.
| // svMP currently uses 2 ms for the M-cell tau_s2 default; | ||
| // Bueno-Orovio et al. (2008), Table 1, reports 4 ms. |
There was a problem hiding this comment.
I suggest rephrasing this as:
| // svMP currently uses 2 ms for the M-cell tau_s2 default; | |
| // Bueno-Orovio et al. (2008), Table 1, reports 4 ms. | |
| // Beware that the default here for the M-cell tau_s2 (X-th entry in | |
| // the list below, 2 ms) is different from the value reported in | |
| // Bueno-Orovio et al. (2008), Table 1, which is 4 ms. |
(replace X-th with 2nd or 3rd depending on which is the right one).
Additionally, perhaps this should be moved to the general class documentation, so that it will be exposed in the Doxygen-generated files. I would suggest placing this in a @warning paragraph.
| * A test-only adapter derived from the concrete model exposes the protected | ||
| * node-local interface needed to run the trajectory experiment. | ||
| * By default, state values are read from CSV columns @c s0 through @c sN. | ||
| * The constructor's state-column mapping can explicitly select other columns | ||
| * when an existing reference file uses a different, documented layout. |
There was a problem hiding this comment.
I feel that this paragraph is a bit hand-wavy, rather than precisely and concretely describe what this test does (it may or may not be AI generated: either way, behave of AI-generated prose, as it very frequently has this issue).
Proposed rephrasing:
* Runs a twitch experiment on the active stress model passed as template
* argument, and compares the simulation result with a reference, trusted
* external simulation loaded from a CSV file. If the simulated result
* doesn't match the reference up to a tolerance, an exception is raised,
* causing the test to fail.
*
* The twitch experiment uses a prescribed calcium transient to trigger
* the contraction in the active stress model (see @ref calcium_at). The
* fiber stretch transient is also prescribed (see @ref fiber_stretch_at)
*
* The active stress model is solved as standalone (i.e. only one system
* of ODEs is solved), not in 3D. Therefore, the test only verifies the
* correctness of the ODE definition and timestepping.
*
* ### Format of the reference solution
*
* The CSV file containing the reference solution is expected to have the
* following format:
* [...insert what the format should be...]
*
* The state variables are assumed by default to have the same order as in
* the concrete model being tested. If necessary, an optional constructor
* argument to this class can define a different ordering.
There was a problem hiding this comment.
Similar considerations apply to this file as to active_stress_test_helpers.h.
In addition to those, it seems to me that the two files aim at the same goal (verification of ODE models against reference solutions), and while they have a similar large-scale structure, there are some small dissimilarities (e.g. the trajectory configuration struct is here but not there, and similar). I suggest removing the dissimilarities as much as possible.
While not technically necessary, this has a few advantages.
- If someone understands the testing infrastructure for one of the two, their understanding transfers immediately to the other (whereas now they could think that it does, but then find out that there are subtle differences, and this may lead them down a wrong path).
- Future migration to a unified setting where
ActiveStressandIonicModelare both under an abstractAuxiliaryStateorODEclass will be easier. As of now it is not trivial to have a unified helper, due to the existing differences betweenActiveStressandIonicModel, and the fact that they do not have a common base class. But I think both of these obstacles will gradually go away.
| # provenance: independent Python Forward-Euler oracle. Equation sources: | ||
| # Nash & Panfilov (2004), doi:10.1016/j.pbiomolbio.2004.01.016, | ||
| # Eq. 22c introduces an excitation-driven active-tension ODE and Eq. 23 | ||
| # supplies its piecewise rate switch. | ||
| # Goktepe & Kuhl (2009), doi:10.1007/s00466-009-0434-z, | ||
| # Eq. 46 reformulates the ODE using dimensional transmembrane potential Phi | ||
| # and resting potential Phi_r; Eq. 47 replaces the piecewise switch with a | ||
| # smooth Gompertz relaxation-rate coefficient. | ||
| # svMultiPhysics-specific adaptation (not documented in either source): | ||
| # substitute intracellular calcium Ca for Phi, with parameter roles | ||
| # eta_T <-> k_sigma, calcium_rest <-> Phi_r, | ||
| # xi_T <-> xi, calcium_crit <-> Phi_bar. | ||
| # The oracle evaluates | ||
| # dT/dt = epsilon(Ca)*(eta_T*(Ca-calcium_rest)-T). | ||
| # Parameters: epsilon_0=0.1, epsilon_i=1.0, xi_T=4e3, eta_T=1e2, | ||
| # calcium_rest=1.25e-4, calcium_crit=8e-4 (slab calibration). | ||
| # Ca transient: double-exponential, c0=1e-4 mM, cmax=9e-4 mM, | ||
| # tau_rise=20ms, tau_decay=50ms, onset=10ms. dt=1ms, 200 steps. | ||
| # Values extracted from the existing test without modification. |
There was a problem hiding this comment.
I appreciate the documentation effort, but as far as I know CSV files do not have a standard comment syntax, meaning that these files might not be readable by other tools. I suggest removing these, and instead placing the information they contain in the documentation for the test case in the associated cpp file.
There was a problem hiding this comment.
I think that the generation of reference solutions should be delegated to trusted third-parties, if possible. Otherwise there's a risk that we are simply comparing two different solvers made by the same people, and both of them might contain the same errors.
I can see that this particular flavor of the Nash-Panfilov model is not the vanilla one (as in the original paper), but it is the slightly modified version from Goktepe and Kuhl, and so the CellML-generated code might not match exactly.
Still, I suggest taking that CellML code (I think it's this one, but do double-check) and applying the modifications there, rather than rewriting the thing from scratch.
A similar suggestion applies to other files as appropriate.
There was a problem hiding this comment.
Also beware that the copyright of CellML-generated code probably still lies with the Cell-ML authors, and so the generated scripts should not contain the header about copyright being held by Stanford University etc.
The CellML model page (at least on the new website) has some licensing and authorship information, and that should be followed in the scripts I think.
There was a problem hiding this comment.
If I understand this right, this file is adapted, if not copied, from this repository. If so, the license statement at the top is not correct (the copyright of at least part of this file is of the original author, not ours).
Having said that, I think for this case it is sufficient to document (e.g. in the README and in the Doxygen documentation) the repository that was used, and the commit hash that was used, rather than including the ready-to-use generator script. I fear that the Python-building-C++-code setup might be quite frail if the system configuration changes (e.g. compilers not found, linking issues, etc.), and it would be an unnecessary hassle to make it more robust.
| # SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the | ||
| # University of California, and others. SPDX-License-Identifier: BSD-3-Clause | ||
|
|
||
| """Regenerate and verify all unit-test reference trajectories.""" |
There was a problem hiding this comment.
What is the verification doing here? I suggest expanding the documentation to explain this more clearly.
Co-authored-by: Michele Bucelli <michelebucelli415@gmail.com>
Resolves #597.
Current situation
The existing regression tests for cellular ODE models use reference solutions generated by svMultiPhysics itself and are generally coupled to larger 3D problems.
This PR adds standalone unit tests for ionic and active-stress models using independently generated reference trajectories.
Release Notes
IonicModelandActiveStress.tests/unitTests/reference_data/.tests/unitTests/reference_generators/.Testing
Code of Conduct & Contributing Guidelines